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We present an exact field theoretical representation of the statistical mechanics of classical hard- 
core Coulomb systems. This approach generalizes the usual sine-Gordon theory valid for point-like 
charges or lattice systems to continuous Coulomb fluids with additional short-range interactions. 
This formalism is applied to derive the equation of state of the restricted primitive model of elec- 
trolytes in the low fugacity regime up to order p 5//2 (p number density). We recover the results 
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I. INTRODUCTION 

The aim of this paper and of the following one is to present a formally exact field theory which allows for the 
calculation of thermodynamic functions of classical hard-core Coulomb systems. It is wclLknown that the grand- 
canonical partition function of the Coulomb gas can be represented by a sine-Gordon actioriJlll. However this mapping 
is applicable only for pointlike charges, or for lattice systemsu If necessary, the short-range repulsion is frequently 
included post-facto by introducing a suitable cutoff in momentum space integrals. In this work, we derive a formally 
exact sine-Gordon held theory whitout the use of any arbitrary cutoff. This off-lattice formalism is used in two 
complementary directions. In-the present-paper we recover the low fugacity expansion of the thermodynamic functions 
obtained years ago by MayeiQ and HagaQ by means of graph resummation techniques. In thej-cpmpanion paper, we 
consider rather the high temperature regime and the results obtained by Stell and Lebowitzt3 in the frame of the 
so-called 7-ordering theory are also recovered. „ 

The simplest sound theory of electrolytes is due to Debye and HuckelEl who showed, nearly eighty years ago, 
that, at least in the low density limit, the potential of mean force ipi2 (f) between two ions of respective charges 
e\ and e2 behaves like eie2 exp(— nr)/r as r — > 00 rather than like the Coulomb potential itself e\eijr, where 
k = (4-7T/3 P^l) 1 ^ 2 is the inverse Debye shielding length. Here /3 = 1/kT (k Boltzmann constant, T temperature), 
Pi is the density of ionic species i and the dielectric constant D of the solvent has been absorbed in the definition of 
the charges. An important consequence of the shielding effect is the non-analyticity of the specific excess osmotic free 
energy f(p) as a function of the mean ionic density p. Actually, in Debye-Huckel theory one obtains f(p) ~ k 3 ~ p 3 / 2 
for p — ► 0. 

The results of Debye and Hiickel are valid only at very low densities and discrepancies between their theory and 
experimental data on real electrolytes have motivated an enormous amount of theoretical works to improve the 
theoretical scheme. A first systematic pertubative. expansion of f(p) in which the density p is taken to be thg 
ordering-,(smai]) parameter wap-proposed by Mayer E Improvements on this seminal work were made later by Haga,Q 
MeeronJaAbeJa and Friedman.ll3 All these works are based on diagrammatic techniques and the more elaborate of 
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them provide a reasonably accurate description of the thermodynamic properties of ionic solutions in the low fugacity 
regime. A monography by-, FriedmanO summarizes the above-mentioned works. 

In a more recent work,l!3 Stell and Lebowitz have proposed a perturbation scheme in which the ordering parameter 
is 7 = (3e 2 where e is the electron charge. Their theory is also based on a sophisticated diagrammatic analysis which 
gives an explicit high temperature expansion of f(p) for symmetric and asymetric electrolytes. 

EdwardsQ seems to have been the first to use the so-called sine-Gordon (SG) transformation in the field of the 
statistical mechanics of classical Coulomb systems as an_alternative to the above-mentioned diagrammatic techniques 
although Kacj3 Siegertjj Hubbardjl3 and StratonovichEJ also pioneered the method in other domains of statistical 
physics or in field theory. In his work, Edwards considers a model of charged hard spheres which allows a clear 
splitting of the pair potential into a long range electrostatic part for which the sine-Gordon transformation applies 
and a repulsive part for which low density virial expansion techniques can be employed. This leads to an intricate 
double expansion in 7 and in p. The sine-Gordan formalism for both classical and quantum Coulomb systems has 
been reviewed recently by Brydges and Martin.Ej 

A decade later, Hubbard and Schofieldlla have shown that a general fluid Hamiltonian with long-range and short- 
range interactions can be mapped onto a reference system with short-range interactions only. Then the cumulapt 
expansion is used to map the original fluid Hamiltonian onto a magneticlike Hamiltonian. Brilliantov and al.cZI 
have explored this route for ionic fluids in order to study the Coulombic criticality. The equation of state ofr-a 
multicomponent system of pointlike ions embedded in a neutralizing background has also been studied by Ortnerlla 
along these lines. Recently Netz and Orlandtll tried to improve on Edwards theory by performing a double 
transform, both on the Coulomb and the hard core parts of the pair potential. As stressed by Brydges and Marti: 
a SG transform of the singular hard core potential and more generally that of a repulsive short range potential such 
that ~ 1/r 12 , is strictly speaking impossible, since these singular potentials do not have a Fourier transform. 

The present work is along the lines of the papers of Brilliantov, Ortner, Netz and Orland. We limit ourselves to 
the case of a symmetric fluid of charged hard spheres with only two species of ions of equal diameters a and carrying 
opposite charges (the so-called restrictive primitive model (RPM) of electrolytesH) . In a first step, we regularize the 
Coulomb potential by a smearing of the chawes-iOver the surface of a sphere of diameter a < a, and therefore give 
a precise meaning to the SG transformation.! 1 5 H 19 I Obviously other kinds of smearing are possible and would lead to 
the same results; in another context an uniform volumic smeacing of the charge has been proposed.! 2 !] This allows us 
to derive rigorously a result which seems to belong to Siegerttl in the general case and which states that the grand- 
partition function of charged hard spheres is equal to the average over a Gaussian measure of the grand-partition 
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function of bare hard spheres in the sine-Gordon field. This is our Eq. (2.1C). In a second step, making a 
connection with liquid theory,K3 we perform a functional expansion of lnHjjg with respect to the sine-Gordon field 



(j) which yields the exact expression (2.14) of the sine-Gordon action £>(</>) of the model. This action involves the 
connected correlation functions of the hard sphere fluids which are supposed to be known. It can be checked that this 
action reduces to the usual sine-Gordon action in the limit of vanishing hard-core diameters. This formalism is very 
handy since it allows to obtain either low fugacity or high temperature expansions of f(p) via cumulant expansions 
for off-lattice Coulombic systems. 

The paper is organized as follows. In Sec. |n|, we derive the generalized sine-Gordon representation of the grand- 
partition function of the RPM model. The low-density expansion of the grand potential (or pressure) of the model 



is obtained in Sec. [HI] up to the order p 5 / 2 . We check that each term of this expansion is actually independent of 



the smearing diameter a. Comparisons with the results of Haga and of Netz-Orland are carried out in Sec. IV. All 
approaches yield identical results at order p 5 / 2 . In addition we explain why the approximate derivation of Netz-Orland 
leads to the correct result at this order but could fail at higher orders in p. Conclusions are drawn in Sec. V. 

II. MODEL AND FORMALISM 
A. The Boltzmann factor 

Throughout this paper we consider the three dimensional (3-D) and symmetric version of the RPM, i.e. a system 
made of AL_ hard spheres of diameter a and charge e and iV_ spheres of the same diameter but with an opposite 
charge — e.El 

With obvious notations, the configurational energy of the model reads as 

P V(r N \r N -) = (t£ (^ i +VHs(r ij )) , (2.1) 
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where e, = ±e and Vf ls (r) denotes the hard core potential. For the moment the ions are supposed to be confined in 
some arbitrary volume l^cM 3 . 

We first note that only configurations (r N + , f N ~ ) of ions without overlaps of the spheres do contribute to the 
canonical (or grand-canonical) partition functions. For these configurations, the charge ±e of any ion of center fl can 
be smeared out uniformly on any spherical surface of diameter < a < a. The interaction energy of two balls of 
charge density r(r) = S(r — a/2)/(ira 2 ), located respectively at point r and r 1 will be noted W T (r — r ') and we have 
obviously 

W T (r-r')= f f d 3 xd 3 yr{\r-x\)-J-^T{\y-r'\) . (2.2) 

Jv Jv f — y\ 

We note that the self-energy £$ = W T (Q)/2 — 1/a of each spherical distribution is a finite quantity for a > 0. 
It follows from the preceding remarks that the electrostatic part of the Boltzmann factor can be written, for any 
configuration (r + , r N ~ ) , as 

ex P E^j = exp(/3iV 7 £ s )exp (-l( n \W T \n)) . (2.3) 

where 7 = fie 2 and N = N+ + N- is the total number of ions. In Eq. (2.3) 

N + N- 

^ = ^^-f,+)-^<5(f-r-_) (2.4) 



n(r) 

i=l 



is the microscopic charge distribution (divided by e) and 

(n\W T \n) = I [d 3 rd 3 r'n(r)Wr{\r-r'\)n(f"). (2.5) 



V JV 



Note that it follows from the positivity of the Fourier transform W T (k) = 47rf(fc) 2 /fc 2 thai the quadratic form (n|W T |n) 
is defin ite p ositive. We can take advantage of this positivity to perform a SG transformllj and reexpress the Boltzmann 
factor ( |2.3| ) as an average over a Gaussian scalar field <f>, i.e. 



N+ N- 



exp (-~(n|W T |n)) =(expyij d 3 x n(x)(f>(x)^ ) Wt = (exp | i ^ <j>(r i+ ) - i (f>(n-) J )w T , (2.6) 



where 4>{x) = ^<j)(x) is a real random field. The precise meaning of the average in Eq. (2J3) is given in Appendix A. 



B. Grand partition function 

For reasons which should become clear below, the grand-canonical ensemble is considerably more handy than the 
canonical one. For simplicity we choose the same chemical potential /1 = [1+ = /i_ for the anions and the cationst3. 
The grand-canonical partition function is given byE2l 

00 00 „N + N- 



{v + ,v_)= E WvWV / ^ + ^-exp(-^(r^,^-)) . (2.7) 

AT . — n AT — n V + / • \ )' JV 
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Ar + =0 Af_=0 

where we have introduced the usual notation v = v± = /3/j, and the activity z = z± = A -3 exp(/3/i) of both species 
(the thermal length A is assumed to be the sam e for the anions and the cations). 
Gathering the intermediate results (2.2), (|2.6|) we get : 



00 00 S JV_ 

■-RPM 



- ) = E JLm^wvL dW ^- exp (-E "Ur v ) ) (2.8) 



N+ N- 



(exp +iJ2^(ri+)-iJ2^(n-) \ )w T ■ (2.9) 



i=l i=l 
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where we have defined renormalized chemical potentials D± = i/± + /3e 2 £s = v + 7 / a an d renormalized activities 
z = z± — exp(P±)/A 3 . 

This last equation can be elegantly rewritten as 



(u + ,u-) = (E + (v + ,u-;i</>, -ut>))w r , 



(2.10) 



where z,+ s _ (i/+,V—;'< 



denotes the grand-canonical partition function of a mixture of two species of equal size 



hard spheres labelled + and — with chemical potentials v + and £L res pectiv ely. The spheres with the label + are in 
the external fiqld icf> whereas those labelled — are in the field —icf>. Eq. ( 2.10 ) is a special case of a more general result 
due to Siegert.0 

In order to get a more explicit expression of the action we perform now a Taylor functional expansion of In S^"^ 



with respect to the activity z± (r) = z e: 
We have, from standard liquid thcorj 



In 



p/a)exp(±i<f>(r)). 
quid theoryHa 

~f s _{v+,v~,V,l3-i4>,-i4>)\ 



s^fL(y,v,v,i3) 



„! a2 



71=1 



Qi ■ ■ -an — i: 



d A \ ■■■d i n 



8z ai {l) ■ ■ •6z an (ri) 



n ( z ^) - z ) 



(2.11) 



i(0= 



In the absence of an external field, we have, for sufficiently large systems zS_ s _{v, v) = S ffs (i/o) , where Vq = i/ + ln2 
(i.e. zq = 2z) and 'E HS \vq) is the grand partition function of a fluid of identical hard spheres. The integral kernels in 
Eq. (2.11) are related to the correlation functions /Iq™^(1, ■■• ,ri) of the hard sphere mixtureil 



8z ai (l) ■ ■■8z cen (n) 



.(i)=* 



In Eq. ( 2.12| ), po = V 1 dla'B HS /dug is the number density of this fluid. Making use of of the identity 

— " = exp (^/ a ) cos ~ 1 



2z 



(2.12) 



(2.13) 



we can rewrite Eq. (2.11) as a functional integral 
S™(i/ +>I /_) 



(exp(-J7M)) Wr =M^ / Vcj> exp(-5[0]) , 



where Mw = J T^<\> exp(— ^ {^| W^T 1 ) is a normalization constant and 



oo 



u[4>] , 



u n [4>] 



n=l 



/ d 3 l ■ ■ ■ d 3 n h (l, ■ ■ ■ , n) [exp (7/0) cos tf>(i) - l] 



i— 1,71 



(2.14) 



(2.15) 



(2.16) 



The above expression of the sine-Gordon like action S[<fi] of the RP M is an exact result; note that S[4>] is an even 
function of the field. In the limit a — > only the term n = 1 of Eq. (2.16) survives and one checks that one recovers 
the usual sine- Gordon action of the Coulomb gas £3 ■— ■ 

A similar result along the Hubbard-Schofield scheme has been used by Brilliantov and alO in their study of 
the criticality of the RPM model. However their approach is developed in Fourier space and without the explicit 
regularization obtained with the smearing of the charges. As a final remark, we note that we might have been 
tempted to perform the Taylor functional expansion, not around z but around z± = zexp(7/a). In that case a 
low-fugacity expansion is valid only at high temperatures since z ~ z requires 7 — > 0. 
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C. The generalised screened potential X T (r) 



The program is now to perform a systematic cwpulant expansion of the expression ( 2.14 ) and to compute the 
cumulants by an extensive use of Wick's theorem. However, as it stands, this expansion will involve cumulants 
which diverge in the thermodynamic limit due to the long range of W T (r) . The same problem arises in low fugacity 
or high temperature-jiiagrammatic expansions of the RPM where one is led to resum classes of diagrams in order 
to get finite results. Bt3E3 Formally it amounts to. introduce a screened (Yukawa) potential. In the field theoretical 
formalism discussed here, a screened (or HartreeEil) field can also be introduced as follows. We denote Uq[(/>] the high 
temperature approximation of U[<f>], i.e. 

U<M = £ / d W) 2 • (2-17) 



2 „ 

Writing now the triviality U[4>] — (U[4>] — Uo[<j>]) + Uo[(j)], we get 

(exp {-U[4>]))w T = ^(exp (- (U[<f>] - U [</>])) , (2.18) 
J\w 

where Ma = J "D<t> exp( — ^ (01^4^7 1 1<^) ) with A = X T ,W T and where X T (r) is a real operator defined by the relation 
Xr(r) -1 = W T (r-) _1 + P07 I (I denotes the identity). 

Using the precise definition of the functional integration given in Appendix A, we obtain 

MX «p(™ / ^i^(l + 7P0 W T (q))) , (2.19) 



M w V 2 J ( 27r ) 



Note that, in Eq. (2.19), we have replaced a scries by an integral, which is valid for large systems. Moreover the 
integral converges for a ^ 0, which is a happy consequence of the regularization of the Coulomb potential via the 
smearing of the charge. 

The Fourier transform X T (q) reads : 

X (a) _ WAq) _ sin 2 ( g q/2) A, 

l + 7Po^(?) (W2) 2 9» + «2=^^ " { } 

with Kq = Attjpq. 

It is easily checked that for a fixed kq and in the limit a — > 0, X T (r) reduces to the familiar screened Yukawa potential, 
i.e. X T (r) ~ exp(— KQr)/r (Vr). Conversely, for a fixed a and in the limit kq — > we have obviously X T (r) = W T (r). 
For arbitrary (a, kq), the large r behavior of X T (r) is determined by the small q behavior of X T (q) ~ 4ir/(q 2 + k,q) 
which implies X T (r) ~ exp(— KQr)/r at large r. The function X T (r) is thus a short range function of r which, however, 
for large Ko, can have a non monotoneous behavior. Expansions of X T (r) at low kq will be given in Sec. |Xj. 

We shall see in the next Section that the use of the screened potential X T (r) instead of the long range (Coulombic) 
potential W T (r) ensures the convergence of the cumulants. 



III. LOW-FUGACITY EXPANSION 



In this section, we use the general result of Eq. (2.14) to perform a systematic low-fugacity expansion of the pressure 
of the RPM model, at the order 5/2 in the density p. 

Although this theory must not depend explicitely of the smearing parameter a, we shall explicitely verify below 
that it is effectively the case for each order of the expansion. 



A. The grand potential 



Using the results of Eqs. ( 2.14 ), ( 2.18 ), ( 2.19 ), and the cumulant theoren£3 we obtain the specific grand potential 
urpm = — 1iiSrpm/V 



urpm{v) = uhsM + 2 



<hT ln(l + 7 po^))-^£ 



(27T) 



1 ^ 



V ^ n\ 

n—1 



(H r ' 



)X T ,c 



(3.1) 
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where H[cj>] =U[4>] — Uq[4>) and (• ■ -)x T ,c denotes a cumulant average. 

Since we are interested in a low-density expansion, we keep only the first two terms U\ and U% of the series ( 2.16 ), 
which is equivalent to take into account all contributions with one and two point correlations functions. Then, 
H[4>] = Ui[4>] + U2[4>] — Uq[4>], and the grand potentiel can be recast in the following form : 



urpm{v) = uhsM + -i + ^2 + O(pl) 



with the following definitions 



d 3 q 

(27T) 3 



In l+7p Wr(g) 



-2 



(H) 



(H) 



V 



2V 



The well-known HS contribution is given by 



2tt 



(3.2) 

(3.3) 
(3.4) 

(3.5) 



uj\ is a generalization, when smearing is taking into acco unt, of the familiar DH expression of the free energy. 
Notice that, unlike the point-like DH approach, the integral (3.3) which defines u)\ is convergent at large k (within 
the conventional theory, an infinite self-energy must be substracted to recover finite results). Since it includes W(q), 
u>i is a function of a (and kq); its expansion in powers of kq is given in Appendix B and reads : 



rf 3 



q 



(2tt)3 



ln( I + ~fp W T (q) 



Pol 
a 



15 



(3.6) 



The first term recast in the form — Kq/127t is reminiscent of the familiar DH contribution to the free energy. 

ll>2 contains averages over the gaussian field <j) which can be obtained from application of Wick's theorem; the 
detailed calculation is reported in Appendix C and one finds 



LJ 2 



^ - ^ + Po [1 - exp(-A /2)] - | [1 - exp(-A /2)] 2 / drtf{r) 
a 2 4 J 

«p(-Ao)f; / drfS {r) 



(3.7) 



^[l-ex P (-A /2)] 2 ^(0) + ^exp(-A ) / ,lrr-Ar) 



Pi 



(2n)l 



where h (0) denotes the 3D Fourier transform of h y '(r). We have introduced in this last expression the dimensionlcss 
potential ip T {r) = jX T (r) and the quantity Ao = tp T (Q) — Expansions of ip T (r) and Ao at low kq are given in 
Appendix B 



27 



A = Vt(0)-^=7 
a 



while 



ip T (r) 



2 7 7 



-«o + Jj- K o a - ^l n l a2 + ( K o) 



r — 7K0 + O(kq) for r<a, 



I I \ 2 ex P(~ K 0r) , 2 ^ 



for r > a 



(3.8) 

(3.9) 
(3.10) 



with go = 2s in h(K oa/2)/ (/coa ). 



Using ( |3 . 8[ ) , ( 3 . 9 ) and (3.10) each contribution to (3.7) can be easily expanded in powers of p$. The first two terms 
give 



PoA 



+ po [1 - exp(-A /2)] 



2„,3 



;Pol 



6 



15 



O(pl) 



(3.11) 



The prefactor of the third term is p\ [1 — exp(— Ao/2)] = ttj 3 Pq + O(p ) and the integral must be splitted into 
two components : 
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r< ff y 15 3 J 

+ 47r 7 2 (a - a) + 8^^7 5/2 pJ /2 (a 2 - ^ 2 ) + O(po) , (3-12) 

drtf(r) = ^^ + 0(l). (3.13) 
Po 

Notice that the last integral is singular in po when po ~~ * 0. Consequently the term under investigation will contribute 
at the order p^ 2 . 

The expansion of the prefactor pgexp(— A ) is straightforward p§exp(— A ) = p\ + 2y / 7r7 3 / 2 ( Og / ' 2 + 0{p%), while 
J r>a drip 2n (r) for n > 2 is related to the exponential integral function E„(z) = dt exp(— zt)t~ n . Indeed one finds 



3C 



-4^ 7 g _Jin_ +OM (3 . 14) 

where the last equality is obtained with the help of the scries representation of the exponential integral functional. 

Gathering all these results, it can be checked that all contributions involving the smearing diameter a exactly cancel 
each other and that we indeed obtain a result independent on a; more precisely we get 

»rpmM = "Po - ^Po 7 V /2 - Po {^f - + f 7 3 + 2^ 3 S( 7 /a)) 



tt0F7 9/2 Po /2 



12 



87 B - 3 + 8 ln(8V^7 1/2 cr) + /2 p^ 2 a - 2n^j 5/2 pl /2 a 2 

- 4 7 rV^7 3/2 Po / V [S(7A0 " T( 7 /o-)] - ^|^ Po /2 Mpo) + G(p 3 ) ■ 
where we have introduced the two following series 

b(7/cr) = > — , 

tri (2ra)!(2n-3) 



(3.15) 



n=2 

\2n+l 



B. Pressure of the RPM model 

In this section we derive the pressure for low density systems at arbitrary temperatures. In order to obtain these 
quantities, we must first performed a transformation from the density p (corresponding to the chemical potential v ) 
to the activity z of the model. For that purpose we use the relation : 

d UHsM _ 47TC7 3 2 3 _ 167TCT 3 

~dzV~ ° " ° ( ° } " " ~3" 

Inserting ( p7| ) in (^15|) we get the expression of lorpm as a function of z and 7 



Po = -^0 T = *o - ^-^ 2 + G(z 3 ) = 2z - ^-^ 2 + 0(z 3 ) . (3.17) 



4V2^7 3/2 ,.3/2 



-w.-wu. t) = -2* - '^3' z^ + 4na 3 (| + (J)* - \ (1)' - 2 S( 7 /<t)) z 2 

+ z 5 / 2 4V2^ 3 / 2 ^-^!(8 7J3 + 81n(8v^ 7 1/2 ^) - 3) + 2 7 7 / 2 a - 2 7 5 / 2 a 2 ) 
, 16A/27r 3 / 2 -7 3 / 2 a 3 , 4\/27r 3 / 2 -7 9 / 2 

+ .5/2 W27T 7 (1 _ 3 p(7/ff) _ T(7/ff)]) _ z5/ 2 ln(2z) lV^J_ + 0(z 3) , (3 _ 18) 
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The density of the system is obtained by the relation p = —zdujRPM{z)/dz which is easily inverted and reads 



Pa V^l 3/2 3/2 . ;; 

z = Y-—2— Po + 7T<T 

v 9/2 



5 7 5/2 (J 2 



+ Po 7 V /2 (~ V (1 ° 7B + 101n ( 8 ^7 1/2 ^) - 7) - 7 V* a 
+ pT^'^ (-§ + 2S( 7 / CT ) + 5T( 7 /a)) - pT Mp^^~ 



o{pI) 



(3.19) 



Using this last equation in (3.18) we obtain the pressure (3Prpm = —Urpm(p,7) 01 the RPM model 



(3Prpm = P 



3/2 



,3/2 , „3 



2tt 



a I — r 7T 

3 \a 



2tt S( 7 /<r) p' 



7rV^7 9/2 + 7s + ln(8V^7 1/2 p 1/2 cr) ) + 37r\/¥7 5/2 cr 2 - 6ttV¥ 7 3/ 2 cr 3 T(7/cr) | /- 
G(p 3 ) 



5/2 



(3.20) 



which can be recasted in the following form : 



ft 2na o k cr o ■ — . 

= p - tj- + ^P 2 + T7^- 2^V S 7 /t - 77^3 
247r 3 1d7t 5127r J 

+ 3 (4 , 7) 3/2^5/2 T(7/a) + 0(p3) 



2 , A 3ft 5 a 3 



(3.21) 



with ft 2 = 47T7ft 2 and ft 2 = 47T7p, the square of the usual inverse DH length. S(-f/a) and Tfa/a) are defined by Eq. 

(B- 



IV. DISCUSSION 



In this section we briefly compare the results obtained above with the classical diagrammatic results of HagaS and 
with the field theoretical approach of Netz-Orlandtil. It will be shown below, as expected, that these three different 
routes yield the same result. 

Haga's expression for the equation of state (Eq. (29) of refQ) is easily compared with Eq. (3.21); both results coincide 



except for the terms proportionnal to ft which differ by a factor of 1/2. In a recent paper 



Bekiranov and Fisheia, 



a slip in Haga's Eq. (25.4) was noted by these authors who pointed out that the last term of Haga's equation should 
read ft 5 a 2 /167r instead of ft 5 a 2 /327r. When this correction is taken into account Haga's results and ours are identical. 

In their papercj, Netz and Orland do not compare explicitely their results with those of Haga. This can bis-done by 
confronting our expression of the grand potential in terms of the fugacity (cf Eq. ( 3.18| )) with Eq. (22) of refEil. Their 
results are given in terms of the hyperbolic sine-integral function Shi(7/cr) and of the incomplete Gamma function 
r(0, 7/er). As shown in Appendix D these functions can be reexpressed in terms of the series S( 7 /it) and T( 7 /<r) 
defined by Eq. ( |3.16 ) and it appears that the two results coincide. Let us now explain, why the approximate theory of 
Netz and Orland gives the correct result at this order. Recall that these authors also perform a Hubbard-Stratonovich 
transformation on the singular hard-core potential. The associate random field is denoted ip( r ) an d the following 
average is needed (Eq. (12) of refl 



(h{l)h(2)} = exp(-w(12)) 



(4.1) 



wher e h(r) = exp(—iip(r) + w(0)/2) and w(r) denotes the hard-core potential. Two consequences result from Eq. 
(4.1). On the one hand, as noted by the authors, exp(— w(12)) (contrary to w) is finite, which regularizes the theory. 
On the other hand, since exp(—w(r)) = 1 — 6{r — a) {9 is the Heaviside function), this procedure amounts to incorpore 
hard-core effects only at the level of the second virial coefficient. It appears (see Eqs. ( C12JC15 ) of Appendix C), that 
this approximation is sufficient for an expansion up to p 5 ^ 2 but will miss some contributions at the next order. 
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V. CONCLUSION 



In conclusion we have proposed in this work a formally exact field theory for hard-core Coulomb systems. This 
approach generalizes the usual sine-Gordon theory valid for pointlikc charges to realistic Coulomb fluids with additional 
short range interactions. Within this formalism we derive the equation of state of the RPM model up to p 5 / 2 . Our 
results confirm the classical diagrammatic expansions of the Mayer-Haga diagrammatic theory. Going to next order is 
perhaps not out of reach although there is a delicate analysis of the relevant contributions to the cumulants to perform. 
Note that an equation of state for pointlike ions as been recently obtained by Ortnenifl using the Hubbard-Schofield 
approach, up to the p 3 contribution. In this latter case, however, the reference system is the ideal gas system which, in 
turn, greatly simplify the calculations. In the companion paper, we derive an equation of state of the RPM in the high 
temperature regime, at any density, by using the formalism developed in the present paper. This two complementary 
limits show the ability of our formulation to tackle in a coherent way the equation of state of Coulomb systems. 

Another problem, which remains a challenge to theory, is the understanding of ionic criticalityO. We believe that 
our formalism might be used to give some insights upon these interesting questions. Work in that direction is currently 
in progress. 
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APPENDIX A: FUNCTIONAL INTEGRATION 



For any operator of the real field 4>(r) we define the average (-4[0])w T by the following relation 



where J T>4> denotes a functional integration and the Gaussian weight P\y T [<ft] is defined as 

exp (-^IW- 1 ^)) 



J T>4> exp i 



\Wr X \<t>)) 



(Al) 



(A2) 



In order to give an unambiguous definition of the measure T><f> and thus a precise meaning to the SG transform we 
henceforth consider a cubic volume V = L with periodic boundary conditions (PBC). The 1/r Coulomb potential 
which enters the configurational energy (2.1) of the RPM must therefore be replaced by the Ewald potentia£3c2l 



exp (iff- r) 



(A3) 



where q = 2-Kn/L (ft = (n x ,n y ,n z ) £ Z 3 ) is a vector of the reciprocal lattice. Recall that E(r) is the periodic 
electrostat ic p qt^utjal of a point charge embedded in a uniform neutralizing background which kills the term q = in 
the series ( A3)B3H This causes a (hopefully) slight difficulty since configurations with N + ^ iV_ are associated with 
the presence of a background which ensures the electric neutrality of the system. The periodical system considered 
here is therefore slightly different from the usual RPM; however, this should make no difference in the thermodynamic 
limit. | — | 
Assuming PBC we thus haveEj 



/. p + oo p+oo 
H J — OO J — CO 



where (^5 and denote respectively the real and imaginary parts of the Fourier component 

d 3 f <fi(r) exp(—iq ■ r) 



(A4) 



(A5) 



of the real field 4>. The infinite product in Eq. (A4) runs over the vectors q ^ of the reciprocal lattice. In fact, 
since, due to the reality of the field (j>, <j>$ = <j>*_^, only half of the vectors has to be considered, for instance those with 



i x > OBB That is what is meant by the subscript ' in Eq. (A4) 
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APPENDIX B: LOW k EXPANSIONS 



In this Appendix, we give the calculation of the expansions of Aq, tp T {f) and u)% with respect to Kq. 



1. Let us begin with A ; we recall that A = Vv(0) — 27/a where ip T (0) = jX T (0). Using Eq. ( p.20| ) we get 



-MO) 



N3^r(g) = — 



sin 2 x , 2 
ax = — 



(2tt) 3 tW.^ ^^si^* 

where the dimensionless parameter £0 = Koa/2 and the function 



sin x 1 



a; 2 +£ 2 l + v{x) 



dx , 



(Bl) 



i/(a;) 



sin 2 a; 



- 1 



(B2) 



satisfies |^(x)| < 1 Vx. Therefore the fraction 1/(1 + ^(.t)) in Eq. (Bl) can be replaced by its series representation 
which yields 



A- T (o) = -£(-ri„(eo) 



n=0 



+oc ■ 2 

sin x 



* 2 + £ 2 



v n (x)dx 



(B3) 



where it should be noted that the integrals X„(£o) in Eq. (B3) are entire functions of £o- Consequently we get 
at order £q 



X T (0) = — [lo(Co) " Zi(fo)] + O«o) , 



where 2o(£o) an d Xi(£o) can be computed by means of the residue theorem 



X i(£o) = 77773 (- 9 + 8 ^o + (12 + 8£o) exp(-2£ ) - (3 + 4£ ) exp(-4£o)) 

1 



16$ 



(B4) 



With the help of Mapple^ one finds 



X r (0) 



7 5 

-no + — nla - 7^K 3 a a 2 + 0(k%) 



(B5) 



from which we deduce Eq. ( p.8[ ). 
2. The same method can be used for the expansion of X T (r). 



d 3 q ~ 1 
Xt ^ = / 7T^3 Xt ^ exp ' f) = — 



-°o -2 



sin xsin(2xr/a) 1 



x(x 2 + Co) 1 + v{x) 



dx 



(B6) 



thus, 



X T {r) = — 

ITT 



+ °° sin 2 xsm(2xr/a) ^ /0 



x(x 2 +£ 2 ) 



(B7) 



is a piecewise function defined for r > a by the expression 
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X<(r) 



2r« 



2r 



2r 



1 - exp( £ ) - exp(-2£ ) sinh(£ — ) + O(^) 



= - - ^ - «o + C(«o) 
a or 



and for r > a by 



X>(r) 



sinh 2 (g ) exp(-^gp) _ smh 2 (K a/2) exp(-re r) 



Ml 



(ko«/2) 2 r 



@ and @ are equivalent to (Q and ( plq ). 
3. We give now the expansion of ui\. We have first 



^ In (l + 7^(5) 



a; 2 In ( 1 + £0-^4- I (/ - r 



which can be written as the sum of three terms 

2 



r 2 n 3 



x z 



?0 



* 2 +£o 2 



sin 2 X 



^a 3 ,/-oo 2:2 + Co 



2„3 51 



(") 



«— 1 /-oo 



7r 2 a J — ' n 

n=2 



x u n (x)dx 



The first integral of the last expression can be easily integrated by parts with the result 

^o 3 



Thus, 



Using 



we obtain 



In ! I ' ^ e ° 



x 2 J x 2 + 



dx 



6 



to 



37T0 3 7T 2 a 3 



1 ^o(eo)-4: 



x v (x)dx 



16Co 



9 + 8£o + (12 + 8£o) exp(-2£ ) - (3 + 4£ ) exp(-4£ )] + 
1 ' 



} , tt£ exp(-2£ ) ^1 , 



1 



2C 3 , ^ a 



3 15 3 
which is the result used in Eq. (EDI). 



4^ ~~ 127T + 240tt ~ 96tt 



APPENDIX C: COMPUTATION OF cj 2 

In this section we give the expansion of loi for low kq . We have 

_ (H) Xr {H 2 ) Xt (H) 2 Xt 



Ll>2 



V 



2V 
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with H[4>] = U t [4>] + U 2 [<i>] - U [<j>], 

Each cumulant implies several averages over the field 4>(r) which can all be deduced from the general average 
X = (cos(Ai</>i) cos(A 2 2 ) cos(A3(/)3))jc t where Ai, A2, A3 are real constants. Starting from 

x = l ( cx P(*(-^i £ i0i + ^e 2 (f>2 + A 3 e 3 ^ 3 )))x T (C2) 



= ^ Y ( ex P {ijd z r(j)(r) (eiAi<5(r- fi) + e 2 A 2 <5(r - r 2 ) + e 3 A 3 5(f - r 3 ))J)x T 



where ej = ±1 and using the fundamental relation for Gaussian integrals 



one gets 

3 



X=~exp(-^^A t 2 ) ^exp(-A 1 A 2 e 1 e 2 X T (12)- A 2 A 3 e 2 e 3 X r (23) - ArA^X,^)) 



cxp ( -^^l^A^ ) (cosh (AiA 2 X T (12)) cosh (A 2 A 3 X T (23)) cosh (AiA 3 X r (13)) 



i=l 

3 



= 1 



/.'' 2 1 a 

(^o 2 )x T - (Uq) 2 Xt _ pi 



2L 3 4 



where A and Vv( r ) are given by Eq. ( p.q ), ( p.9| ) and ( p.lO| ), and we have introduced go(r) = 1 + ^o( r ) 

1 c 



At the lowest order in density the correlation function g^ 1 " 1 (r) = for r < a and 1 otherwise; thus 



Gathering together the two precedent results, we obtain Eq. (3.7). 



(C3) 



( / / </ :! r.MfX(F) )) Xt =cxp(-i / ,/ :! F,</ :! r 2 A' r ln - v, >i( n.)^ ) ) . (C4) 



- sinh(AiA 2 X r (12))sinh(A 2 A 3 X T (23))sinh(AiA3X T (13))) (C5) 

from which all useful formula can be obtained. 

It is convenient to recast the various contribution of lu 2 in the following way 

Po [exp(-A /2) - 1J , (C6) 



dr^ T {r) , (C7) 



{UaUl)x ^ Zph^m^ = I exp(-A /2) / dr^r) , (C8) 

«^ - {U "\J U ^ = -| [1 - 2ex P (-A /2)] MO) (09) 

- exp(-Ao) /" <F [go(r) cosh(VvM) - 1] (CIO) 
= -^[l-exp(-A /2)] 2 / 1 o(0) (Oil) 

- £& exp(-A ) / df 50 (r) [coshty r (r)) - 1] . (C12) 



dr 9 f(r) [cosher)) -1] = ±| drg£\r)tf(r) + f^J dr$\r , ^ 

= \ I - \ L + 5 L dr w • (C13) 
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APPENDIX D: COMPARISON WITH NETZ-ORLAND'S PAPER 



In this Appendix we give some technical details allowing the comparison of Eq. ( 3.1 8| ) with the expression of the 
grand potential in terms of the fugacity calculated by Netz and Orland. 

Eq. (22) of Netz-Orland'paper can be written in our notations (A — > z, a — > a) 



where 



lorpm =2z- h /2 z z ' 2 a z ' 2 - b 2 z 2 a z - b 5/2 z 5 / 2 a 9 / 2 - b ln5/2 z^ 2 a 9 ^ 2 ln(za 3 ) - ■ 



(Dl) 



h/2 — 



-27re 3 - ^ (2e 3 Shi(e) - cosh(e)[4 + 2e 2 ] - 2esinh(e)) 



h/2 = - 



1 50 

2e 3 [r(0, e) + 2 1E + - ln(128^e 3 )] - 2 exp(-e)[2 - e + e 2 ] - — e 3 



V5/2 



(D2) 
(D3) 
(D4) 

(D5) 



with e = 7/(1. 

It can be easily checked that 6 3 / 2 coincides with our coefficient in Eq. ( 3.1SQ . b 2 is recovered using the following 
identity 



e 3 Shi(e) - e 2 cosh(e) - 2 cosh(e) - 2e sinh(e) = 6S(e) - 3e 2 - 2 



(D6) 



where S(e) is the series given by Eq. ( 3.16| ). 

65/2 involves the_^ incomplete Gamma function r(0,e) which is related to exponential integral functions by the 
following relationsEa 



r(0,e)=Ei(e) , 

E„+i(e) = - (exp(-e) - eE„(e)) for n > 1 . 
n 



Thus, from e 3 r(0, e) - exp(-e) (e 2 - e + 2) = -6E 4 (e) we get 
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e 3 r(0, e) - exp(-e) (e 2 - e + 2) = e 3 
where it can be noted that 



lne - 7 B 



G 



2 + 3e-3e 2 + 6^ ( ' 



(m — 3)m! 



E 



(m — 3)to! 



S(e)-T(e). 



(D7) 
(D8) 



(D9) 



(D10) 



Inserting Eqs (D9) and (D10) in Eq. (D4) gives our expression of the z 5 / 2 coefficient in Eq. (3.18). 
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